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The coefficients in van der Poel's equation for calculating the shear modulus of a particulate 

composite have been greatly simplified, making the calculation much less unwieldy. Approximate 
solutions of van der PoeFs equation arc also derived, and it is shown that one of the low order approxi- 
mations is Kerner's equation, or Hashin and Shtrikman's equation for the highest lower hound. The 
kerner approximation is often too low in value when the volume fraction of filler exceeds 0.2, hut it 
can he used to provide further simplification in van der Pod's equation, or it can he used as a first 
approximation in a Newton's method of solution. 
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1. Introduction 

Van der Poel's method [1] l for calculating the 
shear modulus of a particulate eomposite is an excel- 
lent method eapahle of giving accurate results [1, 
2], but because of several disadvantages has not been 
widely used. In the original presentation the calcula- 
tion was complicated. A table of values was therefore 
provided, but this table was limited to materials for 
which Poisson's ratio of the matrix was 0.5. In addition, 
there was an error in the derivation [3], and its effect 
on the results was not known. Recently, however, van 
der Poel's method has been reexamined [4], the error 
corrected, and the method extended for use with matrix 
materials having any value for Poisson's ratio. The 
effect of the error was shown to be minor, and a new 
table of values was provided. Continued studies have 
now resulted in an improved understanding of the 
theory, and have provided significant simplifications 
in the calculation. This work is presented here. 

2. Review of the Method 

A detailed derivation based upon a method de- 
veloped by Frolich and Sack [5] has been given else- 
where [4]. Consider an idealized composite material 
consisting of small spheres imbedded in a matrix. 
The spheres are of approximately the same size, are 
firmly attached to the matrix, and are uniformly dis- 
tributed so that the composite material is macro- 
scopically homogeneous and isotropic. To simplify 



calculations consider a unit length chosen such that 
on the average there is one sphere of filler in each 
spherical volume 47r/3 of unit radius. If a is the radius 
of the filler sphere, then by definition of the unit radius, 
a 3 is equal to the volume fraction <p of the filler. 

Consider a large sphere of homogeneous material. 
Consider also a second sphere of the same size con- 
sisting mostly of the homogeneous material but having 
the following structure in the central region: A filler 
sphere of radius a is located at the origin. Surrounding 
this out to a radius unity is a shell of matrix material, 
which in turn is imbedded in the homogeneous mate- 
rial. The mechanical properties of this homogeneous 
material are assumed to be the same as the average 
macroscopic properties that are sought for the com- 
posite material considered here. 

If the same boundary stresses are applied to each of 
the spheres, it is assumed that the displacements in 
the two spheres will be the same at a distance r^>l, 
except for terms of a high order in 1/r. Thus it is re- 
quired that 

Displacement at r ^> 1 in second sphere 



Displacement at r > 1 in first sphere 



1 



«>3 



Const. 



(i) 



Let the first sphere be subjected at the boundary to 
the deviatoric stress system 



ov = 2rP 2 (Cos 0), T r e = TPi(Cos 0), r r *= 



indicate literature references at the end of this paper. 



= 0. 

(2) 
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where o> is the radial tensile stress, and T r e and T r <t> 
are the tangential shear stresses on the spherical boun- 
ary surface at radius r. The corresponding boundary 
displacements in the r, 0, and (f> directions are 

rp rp 

u r = — P 2 (Cos 0), H0 = — P^Cos 0), u 4) = 0. (3) 

In the above equations T is a stress magnitude, and 
G is the shear or rigidity modulus. P 2 (Cos 6) is the 
second Legendre Polynomial, and P^Cos 6) is its 
derivative with respect to 0. 

The stresses and displacements at any radius r 
within sphere 1 of the homogeneous material are also 
given by eqs (2, 3), but the stresses and displacements 
within sphere 2 having the special structure must be 
found by solving Lame's equations. These solutions 
subject to the requirement of symmetry about the 
6=0. axis, and the requirement that the solutions 
possess the same P 2 symmetry as those of eqs (2, 3), 
have been derived in the previous paper [4] and are 
given by 

u r = [A m r + B m r- 4 + C m r 3 + D m r- 2 ]P 2 (Cos 0) (4) 



us- 



}:A m r-\Bmr-* + 7 -J^C m r* 



+ 



l-2v m 

5 — 4>v m 



D m r-*]p^Coa6) (5) 



o> = 2G„ 



T r6 —G„ 



A m -4,B m r*— =C m ^ 



2(5- y m ) 
5 — 4v m 



D m r~* 



P 2 (Cos 6) (6) 



8 l + 2v m 
A m +-B m r~^ + —^ C m r* 

, 2(l + *m) 3 



5— \v„ 



P'(Cosd) (7) 



where the quantity v is Poisson's ratio. 

In the above solutions the subscript m refers to the 
matrix medium of the region a ^ r ^ 1 in which these 
equations hold. Similar equations holding within the 
filler sphere r ^ a are designated by the subscript /, 
but in order that these solutions be finite at the origin 
the coefficients Bf and Df are set equal to zero. For 
solutions holding in the region r ^ 1 outside the matrix 
shell no subscripts are used. In these solutions the 
coefficients C and D do not appear because of the 
requirement eq (1), and the coefficient A is set equal 
to TjG for consistency with the applied boundary 
stresses, eq (2). 

At the boundaries r=a and r=l. continuity of 
stress and displacement is required. This results in a 
system of eight equations in the seven unknowns 
A/ 9 C/, A m , Bm, C m , D m , B, assuming for the moment 



that the shear modulus G is known. If T is considered 
as an additional unknown, a system of linear homogene- 
ous equations results. In order that the solution of this 
system be nontrivial, or that Af, Cf, A m , etc. can be 
expressed in terms of the quantity 7\ the determinant 
of the matrix of coefficients must vanish. Setting this 
determinant equal to zero evaluates G, the shear 
modulus of the homogeneous material and the quantity 
sought, in terms of the known quantities C/, G m , Vf, v m . 

3. Results 

The determinantal equation for G can be reduced 
to the form 



aX'+pX + y = 0. 



(8) 



where 

a= [4P(7- lOvm) -Sa 7 ] [Q-(8-10v m ) (M - l)a 3 ] 
-126P{M-l)a 3 (l-a 2 ) 2 (9) 

(3 = Z5(l-v m )P[Q-(8-l0v m )(M-l)a*] 

- 15(1 - v m ) [4P(7-lft/ m ) -Sa 7 ] (M-l)a 3 

(10) 

y=-525P(l-^) 2 (M-l)a 3 (11) 

and the quantities P, Q, and S are defined as 

P= (7 + 5^)^ + 4(7-10//) (12) 

Q={8-l(k m )M+7-5v m (13) 

S = 35(7 + 5iy)M(l-i/ m )-P(7 + 5i/ m ) (14) 

M is the ratio Gf/G nu and the quantity X is equal to 
(GIG m ) — 1. Equation (8) has two real roots, one of 
which is negative and extraneous. The positive root 
provides the required value for G. 

Equations (8) through (11), although more compact 
than those given previously [4] are still unwieldy, so 
it becomes worthwhile to search for an approximate 
solution. Such a solution could be obtained as a series 
expansion in powers of a = (p 1 l 3 , but in order to obtain 
sufficient accuracy too many terms would be required. 
A more compact solution is one of the Pade type: 



X = 



MO + All X + Al2 X 2 + 



l + d 1 \ + A 2 X 2 + 

where the parameter X is given by 

,_ (M-l)a 3 _ (M-\)<p 



(15) 



(16) 



To obtain this approximate solution \Q is substi- 
tuted for (M - I) a 3 in eqs (9), (10), (11) for a, f3 and y, 
and these equations are then substituted together with 
eq (15) into eq (8). The unknown constants n u ra 2 , d u 
d 2 , etc. are evaluated by equating to zero the coeffi- 
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cients of X in the resulting expression. Several different 
expansions of the form eq (15) can be obtained, because 
the pairs of constants n 2 , d\\ n&, d 2 ; etc. are interrelated. 
Thus one of the pair members is arbitrary and may be 
set equal to zero. A suitable approximate expression is 

X 4 = 

15(1 -i/, w )\ + 810(1 - v m ) (1 - q 2 ) 2 X 3 

1- (8-10^)\-108(l-a 2 ) 2 [8-10^-3/?/14P]X 3 

(17) 
where R is given by 



R = 4P(7 - 10v m ) -Sa 7 . 
To a first approximation eq (17) becomes 
\5(l-v m )(M-l)<p 
' {8-l0v m )M + 7-5i> m 



(18) 



(19) 



This approximation has been given previously by 
Hashin [6], and by Eshelby [7], who found it by different 
methods. It may be regarded as a generalized form of 
Smallwood's equation [8|. 

To a second approximation eq (17) becomes 



X 2 = 



Q- (8-l0v m )(M-l)<p 



(20) 



It can be shown by algebraic manipulation that this 
approximation is the same as Kerner's equation [9], 
or Hashin and Shtrikman's equation for the highest 
lower bound [TO]. 

Equation (20) provides a means of further simpli- 
fying the exact calculation of X. If the quantities a, /3, 
y are multiplied by the factor 

15(1 -Pm)l[Q- (8-10i/,„) (M-l)a 3 ], 

the following new quantities a , /3' , y' are obtained 

a = 15(1 - v m )R ~ 126P(1 - a 2 ) 2 X 2 (21) 

P' = 525P(1 - v m ) 2 - 15(1 - v m )RXt (22) 

y'=-S2&P(l-v m )*& (23) 

These new quantities a , /3', y' may be substituted 
into eq (8) to calculate X. 

The approximate solution eq (20) can be used with 
Newton's method to obtain a better approximation. 
Let 

F{X) = a'X tl +pX + y' 



a' IX+ 



P 



£ 



2a'/ \4a' 
Then to a first approximation 



or 



F(Xs)=F(X 2 ) + (X N -X 2 )F'(X 2 ) = 0. 
F(Xt) 



X\ — X 2 



F'{X,) 



(24) 



(25) 



(26) 



Evaluating the quantities F{Xi) and F' (X-z) gives 



X N = 



x 2 + 



U6P(l-a 2 ) 2 Xi 



525P(l-v m ) 2 +l5(l-p m )RX 2 -252P(l-a 2 ) 2 X 2 l 

(27) 

The function F(X) is a parabola symmetric about 
the axis X = — /3'/2a f . At the point where F(X) inter- 
sects this axis F(X) is negative. As X increases 
positively F(X) becomes less negative and is zero 
where it intersects the X axis. The value of F(X 2 ) is 



F(X 2 ) = -l26P(l-a 2 ) 2 X» 



(28) 



and is negative, so that eq (20) provides a value of G 
less than the true value. On the other hand the value 
ofF(X N ) is 



F(X N )=a'(X N -X 2 ) 2 



(29) 



and is positive. If X 2 is a sufficiently good approxima- 
tion, Xn will be an improved approximation, slightly 
exceeding the true value. 

4. Discussion 

The exact solution of eq (8) and the approximate 
solutions eqs (17) and (27) have been studied for 
volume fractions of filler <p between 0. and 1.0, using 
the following values for the parameters Vf, v m , and M\ 



Vf= 0.1; v 
Vf = 0.1; v 
*>/=0.25; v 
v s = 0.25; v 
i// =0.25; v 
v f =0.35; v 
Vf= 0.35; v 



= 0.4; M = 30, 100 

= 0.5; M = 1000 

= 0.3; M=10 

= 0.4; M = 30, 100 

= 0.5; M = 100, 10\ 10 4 , 7 X 10 4 , 

= 0.4; M = 30, 100 

= 0.5; M = 1000. 
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Figure 1. Plot of relative shear modulus G/G m versus volume 
fraction of filler for the case G f /G m = 30, ^ f = 0.25, ^ m = 0.40. 

Solid curve, exact solution of eq (8); Dash curve 2, X t approximation (Kerner's equation), 
eq (20); Dash curve 3, X :i approximation, eq (17); Dash curve 4, X* approximation, eq (17); 
Dash curve N, Newton's method approximation X\. eq (27). 
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FIGURE 2. f7o£ o/ relative shear modulus G/G m versus volume 
fraction of filler for the case G f /G m = 70,000, ^ f =0.25, ^ m = 0.50. 

Solid curve, exact solution of eq (8); Dash curve 2, X-> approximation (Kerner's equation), 
eq (20); Dash curve 3, Xa approximation, eq (17); Dash curve 4, X» approximation, eq (17); 
Dash curve N, Newton's method approximation X\, eq (27). 



Plots selected from these calculations are given for 
two representative composite materials in figures 
1 and 2. In figure 1 the dependence of relative shear 
modulus GIG m on the volume fraction of filler <p is 
shown, as calculated from eqs (8), (17) and (27). The 
parameter values used were Vf = 0.25, v m = 0.4 and 
M = 30, and are representative of a composite consist- 
ing of small glass spheres imbedded in a rigid epoxy 
matrix. For this example it is seen that the plots for the 
X 4 and X.x approximations lie reasonably close to either 
side of the plot for the exact value, but that the plotted 
values for the X 2 approximation are too low when the 
volume fraction <p exceeds 0.4. 

Figure 2 is a similar plot of values calculated using 
the parameter values z^=0.25, ^ m = 0.50 and 
M =70,000. These values are representative of a 
composite consisting of small glass spheres imbedded 
in a matrix of lightly vulcanized rubber. In this ex- 
ample the Xn approximation calculated from eq (27) 
gives values close to the exact values calculated from 
eq (8), when <p is less than 0.4. The approximate values 
X 2 , Xs, Xa are good only for values of <p up to 0.2. 

Equation (16) for the parameter X may be expressed 
in the following form 



When i> m = 0.4, M = 30, <p = 0.4, the value of X is 
k = 0.09, and when v m = 0.5, M = 70,000, <p = 0.4, the 
value of X is A = 0.13. Thus for values of ip greater 
than 0.4 the terms in eq (17) representing the suc- 
cessive approximations X 2 , X 3 , X 4 are relatively large, 
indicating that this series solution does not converge 
very rapidly, as evidenced in figures 1 and 2. The 
approximations provided by eq (17) are probably 
reliable only for values of <p up to 0.2. In this connec- 
tion it is appropriate to recall that the approximation 
X 2 is equivalent to the widely used equation of Kerner 
[9], or Hashin and Shtrikman's equation for the highest 
lower bound [10]. 

It has been pointed out previously that values of 
G/G,,, calculated from van der Poel's equation are 
insensitive to values of vj. The extent of this insensi- 
tivity can be seen from eq (17). The values of R and P 
are moderately dependent on the values of i>/, and these 
quantities first appear in the coefficient of X 3 in the 
denominator of eq (17). In other words the approxima- 
tion Xa is the first one to be affected by the value of Vf. 

Although in the examples of figures 1 and 2 the 
approximation of eq (27) gives reasonably accurate 
values, there are situations when eq (27) fails because 
the approximation X 2 is not good enough. This occurs 
when for a sufficiently high value of <p the coefficient 
/3', eq (22), changes from a positive to a negative 
value. The axis of the parabola F(X) = 0. then inter- 
sects the X axis at the positive point X = (— fi')l2a . 
This point will usually be near X 2 , and as <p increases 
will exceed it. In either event the slope F' (X 2 ) will 
be too low or even negative, and the calculated value 
of Xn will be grossly inaccurate. 

In the studies using the parameter values listed 
above, it was found that a was always positive and 
y' was always negative, but when <p exceeded 0.5, 
p' changed sign and became negative. For values of 
(p exceeding 0.5 the values of Xn were reasonably 
accurate when M was less than 100, but when M ex- 
ceeded 100 the values of Xn were grossly inaccurate 
or negative. 

It is fortunate that eq (27) gives a good approxima- 
tion for values of (p up to 0.4 or 0.5. At values of <p 
greater than 0.5 the assumptions upon which the theory 
is based are no longer valid. For instance, the packing 
fraction of spheres of equal size arranged on a cubic 
lattice is 0.52, and the packing fraction of equal size 
spheres in hexagonal close packed arrangement is 
0.74. As a practical matter it is difficult to prepare a 
particulate composite material with a volume fraction 
of filler (p exceeding 0.5 without introducing excessive 
void content. However Schwarzl [2] and van der Poel 
[1] report that the theory agrees well with experimental 
data for values of (p up to 0.5. It seems then that the 
calculations of the theory have an empirical validity 
beyond the theoretical expectations. 
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